/* Copyright(C) 1999, Swedish Institute of Computer Science */

#include "fd.h"
#include "dvars.h"

#if MULTI_SP_AWARE
#define add_ddd(A1,A2,A3) add_ddd(HIDDEN_ARG, A1,A2,A3)
#define sub_ddd(A1,A2,A3) sub_ddd(HIDDEN_ARG, A1,A2,A3)
#define div_dud(A1,A2,A3) div_dud(HIDDEN_ARG, A1,A2,A3)
#define div_dll(A1,A2) div_dll(HIDDEN_ARG, A1,A2)
#define div_dld(A1,A2,A3) div_dld(HIDDEN_ARG, A1,A2,A3)
#define div_dldf(A1,A2,A3) div_dldf(HIDDEN_ARG, A1,A2,A3)
#define div_dldc(A1,A2,A3) div_dldc(HIDDEN_ARG, A1,A2,A3)
#define scalar_product_init(A1,A2) scalar_product_init(HIDDEN_ARG, A1,A2)
#define scalar_product_le(A1,A2) scalar_product_le(HIDDEN_ARG, A1,A2)
#define scalar_product_ge(A1,A2) scalar_product_ge(HIDDEN_ARG, A1,A2)
#endif /* MULTI_SP_AWARE */

/* operations on double word entities, for scalar_product */

/* Some "numbers" and their representation:

   0          : < 0, 0, 0>
   0-sup      : <-1, 0, 0>
   0+sup      : < 1, 0, 0>
  -1          : <-1,-1,-1>
  -1-sup      : <-2,-1,-1>
  -1+sup      : < 0,-1,-1>
   1          : < 0, 0, 1>
   1-sup      : <-1, 0, 1>
   1+sup      : < 1, 0, 1>

 */

struct dword {
  SP_WORD infw;
  union {
    SP_UWORD asuword;
    SP_WORD asword;
  } msw;
  union {
    SP_UWORD asuword;
    SP_WORD asword;
  } lsw;
};

typedef struct dword *Dword;

#define issmall_d(x)					\
    (((x)->infw==0 &&					\
      (x)->msw.asword==0 &&				\
      (x)->lsw.asword>=0 && (x)->lsw.asword<HighInt) ||	\
     ((x)->infw==-1 &&					\
      (x)->msw.asword==-1 &&				\
      (x)->lsw.asword<0 && (x)->lsw.asword>=-HighInt))	\

#define isfinite_d(x)				\
    (((x)->infw==0 &&				\
      (x)->msw.asword>=0) ||			\
     ((x)->infw==-1 &&				\
      (x)->msw.asword<0))			\

#define islong_d(x)				\
    (((x)->infw==0 &&				\
      (x)->msw.asword==0 &&			\
      (x)->lsw.asword>=0) ||			\
     ((x)->infw==-1 &&				\
      (x)->msw.asword==-1 &&			\
      (x)->lsw.asword<0))			\

#define cmp_deqd(X,Y)				\
    (((X)->infw==(Y)->infw) &&			\
     ((X)->lsw.asword==(Y)->lsw.asword) &&	\
     ((X)->msw.asword==(Y)->msw.asword))	\

#define cmp_dltd(x,y)				\
    (((x)->infw==(y)->infw) ?			\
      (((x)->msw.asuword==(y)->msw.asuword) ?	\
       ((x)->lsw.asuword < (y)->lsw.asuword) :	\
       ((x)->msw.asuword < (y)->msw.asuword)) :	\
     ((x)->infw < (y)->infw))			\

#define cmp_dgtd(x,y)				\
    (((x)->infw==(y)->infw) ?			\
      (((x)->msw.asuword==(y)->msw.asuword) ?	\
       ((x)->lsw.asuword > (y)->lsw.asuword) :	\
       ((x)->msw.asuword > (y)->msw.asuword)) :	\
     ((x)->infw > (y)->infw))			\

#define cmp_deqz(x)				\
    ((x)->infw==0 &&				\
     (x)->msw.asword==0 &&			\
     (x)->lsw.asword==0)			\

#define cmp_dltz(x)				\
    ((x)->infw < 0)				\

#define cmp_dgtz(x)				\
    (((x)->infw==0) ?				\
      (((x)->msw.asuword==0) ?			\
       ((x)->lsw.asuword > 0) :			\
       ((x)->msw.asuword > 0)) :		\
     ((x)->infw > 0))				\

#define cmp_dnez(x) (!cmp_deqz(x))

#define cmp_dgez(x) (!cmp_dltz(x))

#define cmp_dlez(x) (!cmp_dgtz(x))

static void 
add_ztd(TAGGED x,Dword y)
{
  if (TagIsSmall(x)) {
    y->lsw.asword = GetSmall(x);
    y->msw.asword = (Tltz(x) ? -1 : 0);
  } else {
    TAGGED *s = TagToSTR(x);
    y->lsw.asuword = s[1];
    if (s[0]==BIGNUM_HEADER)
      y->msw.asword = ((SP_WORD)s[1]<0 ? -1 : 0);
    else
      y->msw.asuword = s[2];	/* a two-word bignum */
  }
  y->infw = (y->msw.asword<0 ? -1 : 0);
}

/* carry logic ripped off from bn_add() */
static void 
add_ddd MAGIC (HIDDEN_PROTO Dword x,Dword y,Dword z)
{
  struct dword sum_mem;
  Dword sum = &sum_mem;
  SP_UWORD carry;

  sum->lsw.asuword = x->lsw.asuword + y->lsw.asuword;
  carry = (sum->lsw.asuword < x->lsw.asuword);
  sum->msw.asuword = x->msw.asuword + y->msw.asuword + carry;
  if (carry) {
    carry = (sum->msw.asuword <= x->msw.asuword);
  } else {
    carry = (sum->msw.asuword <  x->msw.asuword);
  }
  sum->infw = x->infw + y->infw + carry;
  if ((x->infw<0)==(y->infw<0) && (y->infw<0)!=(sum->infw<0))
    fd.fd_overflow = TRUE;
  *z = *sum;
}

/* carry logic ripped off from bn_subtract() */
static void 
sub_ddd MAGIC (HIDDEN_PROTO Dword x,Dword y,Dword z)
{
  struct dword difference_mem;
  Dword difference = &difference_mem;
  SP_UWORD carry;

  difference->lsw.asuword = x->lsw.asuword - y->lsw.asuword;
  carry = (difference->lsw.asuword > x->lsw.asuword);
  difference->msw.asuword = x->msw.asuword - y->msw.asuword - carry;
  if (carry) {
    carry = (difference->msw.asuword >= x->msw.asuword);
  } else {
    carry = (difference->msw.asuword >  x->msw.asuword);
  }
  difference->infw = x->infw - y->infw - carry;
  if ((x->infw<0)!=(y->infw<0) && (y->infw<0)==(difference->infw<0))
    fd.fd_overflow = TRUE;
  *z = *difference;
}

static void 
sub_dld(Dword x,long y,Dword z)
{
  struct dword difference_mem;
  Dword difference = &difference_mem;
  SP_UWORD carry;
  SP_UWORD sign_extension = (y<0 ? -1 : 0);

  difference->lsw.asuword = x->lsw.asuword - (SP_UWORD)y;
  carry = (difference->lsw.asuword > x->lsw.asuword);
  difference->msw.asuword = x->msw.asuword - sign_extension - carry;
  if (carry) {
    carry = (difference->msw.asuword >= x->msw.asuword);
  } else {
    carry = (difference->msw.asuword >  x->msw.asuword);
  }
  difference->infw = x->infw - sign_extension - carry;
  *z = *difference;
}

static void 
sub_zdd(Dword y,Dword z)
{
  SP_UWORD carry;

  z->lsw.asuword = - y->lsw.asuword;
  carry = (z->lsw.asuword > 0);
  z->msw.asuword = - y->msw.asuword - carry;
  if (carry) {
    carry = 1;
  } else {
    carry = (z->msw.asuword >  0);
  }
  z->infw = - y->infw - carry;
}

static void 
add_dud(Dword x,SP_UWORD y,Dword z)
{
  struct dword sum_mem;
  Dword sum = &sum_mem;
  SP_UWORD carry;

  sum->lsw.asuword = x->lsw.asuword + y;
  carry = (sum->lsw.asuword < x->lsw.asuword);
  sum->msw.asuword = x->msw.asuword + carry;
  if (carry) {
    carry = (sum->msw.asuword <= x->msw.asuword);
  } else {
    carry = (sum->msw.asuword <  x->msw.asuword);
  }
  sum->infw = x->infw + carry;
  *z = *sum;
}

static void 
sub_dud(Dword x,SP_UWORD y,Dword z)
{
  struct dword difference_mem;
  Dword difference = &difference_mem;
  SP_UWORD carry;

  difference->lsw.asuword = x->lsw.asuword - y;
  carry = (difference->lsw.asuword > x->lsw.asuword);
  difference->msw.asuword = x->msw.asuword - carry;
  if (carry) {
    carry = (difference->msw.asuword >= x->msw.asuword);
  } else {
    carry = (difference->msw.asuword >  x->msw.asuword);
  }
  difference->infw = x->infw - carry;
  *z = *difference;
}


static void 
mul_uud(SP_UWORD x,SP_UWORD y,Dword z)
{
  int carry = 0;
  SP_UWORD lsw;
  SP_UWORD x0 = LOWER_UHALF(x);
  SP_UWORD x1 = (UHALF)(x>>HALFSHIFT);
  SP_UWORD y0 = LOWER_UHALF(y);
  SP_UWORD y1 = (UHALF)(y>>HALFSHIFT);
  SP_UWORD p01 = (SP_UWORD)x0*(SP_UWORD)y1;
  SP_UWORD p10 = (SP_UWORD)x1*(SP_UWORD)y0;
  
  z->lsw.asuword = (SP_UWORD)x0*(SP_UWORD)y0;
  z->msw.asuword = (SP_UWORD)x1*(SP_UWORD)y1;
  lsw = z->lsw.asuword + (p01<<HALFSHIFT);
  if (lsw<z->lsw.asuword)
    carry++;
  z->lsw.asuword = lsw;
  lsw = z->lsw.asuword + (p10<<HALFSHIFT);
  if (lsw<z->lsw.asuword)
    carry++;
  z->lsw.asuword = lsw;
  z->msw.asuword += (p01>>HALFSHIFT) + (p10>>HALFSHIFT) + carry;
  z->infw = 0;
}

static void 
mul_ltd(SP_WORD x,TAGGED ty,Dword z)
{
  z->infw = 0;
  z->msw.asword = 0;
  z->lsw.asword = 0;
  if (ty==Inf) {
    z->infw = (x>0 ? -1 : 1);	/* N.B. just +- infinity */
  } else if (ty==Sup) {
    z->infw = (x>0 ? 1 : -1);	/* N.B. just +- infinity */
  } else {
    SP_WORD y = GetSmall(ty);
  
    if (x>>HALFSHIFT==x>>(HALFSHIFT-1) &&
	y>>HALFSHIFT==y>>(HALFSHIFT-1)) {
      z->lsw.asword = x*y;
      z->msw.asword = (z->lsw.asword<0 ? -1 : 0);
      z->infw = (z->lsw.asword<0 ? -1 : 0);
    } else {
      int negated = 0;
      
      if (x<0)
	x = -x, negated ^= 1;
      if (y<0)
	y = -y, negated ^= 1;

      mul_uud(x,y,z);
      if (negated)
	sub_zdd(z,z);
    }
  }
}

/* Roll-your-own division and remainder algorithms.
   Assume everything is unsigned.

   Consider X/Y, X = (A<<32)+B, |X|<=63, |Y|<=25.

   Let:    Q1 = (1<<32)/Y       |Q1|<=31
           R1 = (1<<32)%Y       |R1|<=25
	   Q2 = B/Y             |Q2|<=6
	   R2 = B%Y             |R2|<=25

   Case 1 (A=0): X/Y = Q2, X%Y = R2.

   Case 2 (A>0): X/Y = A*Q1 + Q2 + (A*R1 + R2)/Y
                 X%Y = (A*R1 + R2)%Y
*/

/* x is nonnegative! */
static void 
div_dud MAGIC(HIDDEN_PROTO Dword x,SP_UWORD y,Dword z)
{
  struct dword dividend = *x;
  struct dword term;
  SP_UWORD q1, r1, q2, r2;

  z->infw = x->infw;
  if (y==1) {
    z->lsw.asuword = x->lsw.asuword;
    z->msw.asuword = x->msw.asuword;
  } else {
    z->lsw.asuword = 0;
    z->msw.asuword = 0;
    /* note: dividend is two times too small */
    q1 = (((SP_UWORD)((-1L)<<(WORDSIZE-1)))/y)<<1;
    r1 = (((SP_UWORD)((-1L)<<(WORDSIZE-1)))%y)<<1;
    if (r1>=y) {
      q1++;
      r1 -= y;
    }
    q2 = dividend.lsw.asuword/y;
    r2 = dividend.lsw.asuword%y;
    while (dividend.msw.asuword>0) {
      mul_uud(q1,dividend.msw.asuword,&term);
      add_ddd(z,&term,z);
      add_dud(z,q2,z);
      mul_uud(r1,dividend.msw.asuword,&dividend);
      add_dud(&dividend,r2,&dividend);
      q2 = dividend.lsw.asuword/y;
      r2 = dividend.lsw.asuword%y;
    }
    add_dud(z,q2,z);
  }
}

/* round towards zero */
static void 
div_dld MAGIC(HIDDEN_PROTO Dword x,SP_WORD y,Dword z)
{
  if (y==1) {
    *z = *x;
  } else {
    int negx=0, negy=0;
    if (cmp_dltz(x)) {
      negx=1;
      sub_zdd(x,x);
    }
    if (y<0) {
      negy=1;
      y = -y;
    }
    div_dud(x,y,z);
    if (negx)
      sub_zdd(x,x);
    if (negx^negy)
      sub_zdd(z,z);
  }
}

static SP_WORD 
div_dll MAGIC(HIDDEN_PROTO Dword x,SP_WORD y)
{
  struct dword z_mem;
  Dword z = &z_mem;

  div_dld(x,y,z);
  return z->lsw.asword;
}

/* round towards -INF (floor) */
static void 
div_dldf MAGIC(HIDDEN_PROTO Dword x,SP_WORD y,Dword z)
{
  if (y==1) {
    *z = *x;
  } else {
    int negx=0, negy=0;
    if (cmp_dltz(x)) {
      negx=1;
      sub_zdd(x,x);
    }
    if (y<0) {
      negy=1;
      y = -y;
    }
    if (negx^negy)
      sub_dld(x,1-y,x);
    div_dud(x,y,z);
    if (negx^negy)
      sub_dld(x,y-1,x);
    if (negx)
      sub_zdd(x,x);
    if (negx^negy)
      sub_zdd(z,z);
  }
}

/* round towards +INF (ceiling) */
static void 
div_dldc MAGIC(HIDDEN_PROTO Dword x,SP_WORD y,Dword z)
{
  if (y==1) {
    *z = *x;
  } else {
    int negx=0, negy=0;
    if (cmp_dltz(x)) {
      negx=1;
      sub_zdd(x,x);
    }
    if (y<0) {
      negy=1;
      y = -y;
    }
    if (negx==negy)
      sub_dld(x,1-y,x);
    div_dud(x,y,z);
    if (negx==negy)
      sub_dld(x,y-1,x);
    if (negx)
      sub_zdd(x,x);
    if (negx^negy)
      sub_zdd(z,z);
  }
}

/* x is nonnegative! */
static SP_UWORD 
mod_duu_internal(Dword x,SP_UWORD y)
{
  struct dword dividend = *x;
  SP_UWORD r1, r2;

  if (y==1)
    return 0;
  /* note: dividend is two times too small */
  r1 = (((SP_UWORD)((-1L)<<(WORDSIZE-1)))%y)<<1;
  if (r1>=y) {
    r1 -= y;
  }
  r2 = dividend.lsw.asuword%y;
  while (dividend.msw.asuword>0) {
    mul_uud(r1,dividend.msw.asuword,&dividend);
    add_dud(&dividend,r2,&dividend);
    r2 = dividend.lsw.asuword%y;
  }
  return r2;
}

/* see adjust_bounds_gcd for x<0 case! */
static SP_UWORD 
mod_duu(Dword x,SP_UWORD y)
{
  SP_UWORD m;

  if (cmp_dgez(x))
    return mod_duu_internal(x,y);
  sub_zdd(x,x);
  m = mod_duu_internal(x,y);
  sub_zdd(x,x);
  return (m==0 ? 0 : y-m);	/* modulo, not remainder */
}

typedef long TERM;

struct linear_data {
  void (SPCDECL *destructor)(void *);
  void (SPCDECL *daemon)(HIDDEN_PROTO
			 Argtype,void *,SP_globref,TAGGED); /* (w,handle,attr_ref,global) */
#if MULTI_SP_AWARE
  SPEnv *spenv;
#endif /* MULTI_SP_AWARE */

  SP_globref refbase;
  int op;			/* 1(#>=), 2(#>=), 3(#=), 4(#\=) */
  long stamp;			/* increases up to backtracking */
  int units;			/* maintained incrementally */
  int nonunits;			/* maintained incrementally */
  long gcdall;			/* valid during GCD rules */
  int prunings;			/* counter */
  int nvars;			/* #terms */
  int ntargets;			/* #terms that may be targets */
  int overflow;			/* set by the daemon */
  struct dword rhs_mem;
  struct dword bige_mem;
  struct dword bigf_mem;
  Dword rhs;
  Dword bige;
  Dword bigf;
  TERM *heap;
  int *vheap;
  int heapsize;
  TERM *target;
  Dvar dvar;
  struct {
    Dword cmin;		/* min(ai*xi) */
    Dword cmax;		/* max(ai*xi) */
    Dword interval;	/* cmax-cmin */
    long *coeff;		/* ai */
    long *abscoeff;		/* |ai| */
    long *gcd;			/* valid during GCD rules */
  } term;
};


#define SV(I)   (pdata->target[I])
#define COEFF(t) (pdata->term.coeff[t])
#define ABSCOEFF(t) (pdata->term.abscoeff[t])
#define GCD(t) (pdata->term.gcd[t])
#define DVAR(t) (pdata->dvar+(t))
#define RefAttr(T) (pdata->refbase + 2*(T))
#define RefVar(T) (pdata->refbase + 2*(T) + 1)
#define CMIN(t) (pdata->term.cmin+(t))
#define CMAX(t) (pdata->term.cmax+(t))
#define INTERVAL(t) (pdata->term.interval+(t))

static void SPCDECL linear_destructor(void *pdata_v)
{
  struct linear_data *pdata = (struct linear_data *)pdata_v;
  FD_SETUP_SPENV(pdata->spenv)

  SP_free_globrefs(pdata->refbase,pdata->nvars<<1);
  SP_free(pdata);
}

#define SWAP(I,J)				\
{						\
  TERM vi = pdata->heap[I];			\
  TERM vj = pdata->heap[J];			\
  pdata->heap[I] = vj;				\
  pdata->heap[J] = vi;				\
  pdata->vheap[vi] = (J);			\
  pdata->vheap[vj] = (I);			\
}

static void 
spheapify(struct linear_data *pdata,
	  int i)
{
  TERM *heap = pdata->heap;
  Dword key = pdata->term.interval;
  
  for (;;) {
    int l = (i<<1)+1;
    int topmost = i;
    if (l<pdata->heapsize && cmp_dgtd((key+heap[l]),(key+heap[topmost])))
      topmost = l;
    if (l+1<pdata->heapsize && cmp_dgtd((key+heap[l+1]),(key+heap[topmost])))
      topmost = l+1;
    if (topmost==i)
      break;
    SWAP(i,topmost);
    i = topmost;
  }
}

static INLINE long gcd (long c1,long c2)
{
  if (c1 > 1 && c2 > 1) {
    while (c1 != 0 && c2 != 0) {
      c1 %= c2;
      if (c1 != 0)
	c2 %= c1;
    }
    return c1+c2;
  } else {
    return 1;
  }
}

static void 
refresh_gcd(struct linear_data *pdata)
{
  int ntargets = pdata->ntargets;
  int units = pdata->units;
  int nonunits = pdata->nonunits;
  long g;
  int i;

  if (nonunits==0)
    pdata->gcdall = 1;
  else if (units>1) {
    pdata->gcdall = 1;
    for (i=ntargets-1; i>=0; i--)
      GCD(SV(i)) = 1;
  } else if (units==1) { /* GCD(i)=gcdall=1 except for one i */
    int ix=0;			/* avoid false alarm */
    pdata->gcdall = 1;
    g = -1;
    for (i=0; i<ntargets; i++) {
      TERM elt = SV(i);
      long ac = ABSCOEFF(elt);

      GCD(elt) = 1;
      if (cmp_deqz(INTERVAL(elt)))
	;
      else if (ac==1)
	ix = elt;
      else if (g== -1)
	g = ac;
      else
	g = gcd(g,ac);
    }
    GCD(ix) = g;
  } else {	/* for each i, compute GCD of all coefficients except i */
    int i0;
    for (i0=0; i0<ntargets && cmp_deqz(INTERVAL(SV(i0))); i0++)
      ;
    g = ABSCOEFF(SV(i0));
    for (i=i0+1; i<ntargets; i++) {
      TERM elt = SV(i);
      
      if (cmp_dgtz(INTERVAL(elt))) {
	GCD(elt) = g;
	g = gcd(g,ABSCOEFF(elt));
      }
    }
    pdata->gcdall = g;
				/* now, SV(i)->gcd = gcd(a[0] ... a[i-1]) */
    for (i=ntargets-1; i>=0 && cmp_deqz(INTERVAL(SV(i))); i--)
      ;
    g = ABSCOEFF(SV(i));
    for (i--; i>i0; i--) {
      TERM elt = SV(i);
      
      if (cmp_dgtz(INTERVAL(elt))) {
	GCD(elt) = gcd(g,GCD(elt));
	g = gcd(g,ABSCOEFF(elt));
      }
    }
    GCD(SV(i0)) = g;
				/* now, SV(i)->gcd = gcd({a[j] | i!=j}) */
  }
}

static void SPCDECL 
linear_daemon MAGIC (HIDDEN_PROTO
		     Argdecl,
		     void *vdata,
		     SP_globref attr_ref,
		     TAGGED global)
{
  struct linear_data *pdata = (struct linear_data *)vdata;
  TAGGED tstate;
  int ar, state_stamp;

  tstate = RefMutable(CTagToArg(global,1));
  ar = Arity(TagToHeadfunctor(tstate));
  state_stamp = GetSmall(CTagToArg(tstate,ar));
  if (pdata->stamp==state_stamp) { /* incremental */
    int j = (attr_ref - pdata->refbase)>>1;
    long c = COEFF(j);
    long ac = ABSCOEFF(j);
    Dvar dv = DVAR(j);
    struct dword cminj_mem;
    Dword cminj = &cminj_mem;
    struct dword cmaxj_mem;
    Dword cmaxj = &cmaxj_mem;
    TERM elt;

    dvar_init(dv, attr_ref, attr_ref+1);
    if (c>0) {
      mul_ltd(c,dvar_min_t(dv),cminj);
      mul_ltd(c,dvar_max_t(dv),cmaxj);
    } else {
      mul_ltd(c,dvar_max_t(dv),cminj);
      mul_ltd(c,dvar_min_t(dv),cmaxj);
    }
    if (cmp_deqd(cminj,CMIN(j)) && cmp_deqd(cmaxj,CMAX(j)))
      return;			/* self-invocation */
    (void)fd_daemon_copy_state(w,&global);
    pdata->stamp++;
    add_ddd(pdata->bigf,CMIN(j),pdata->bigf);
    sub_ddd(pdata->bige,CMAX(j),pdata->bige);
    *CMIN(j) = *cminj;
    *CMAX(j) = *cmaxj;
    sub_ddd(cmaxj,cminj,INTERVAL(j));
    sub_ddd(pdata->bigf,cminj,pdata->bigf);
    add_ddd(pdata->bige,cmaxj,pdata->bige);
    if (dvar_is_integer(dv)) {
      if (ac==1)
	pdata->units--;
      else
	pdata->nonunits--;
    }
    if (fd.fd_overflow) {
      pdata->overflow = TRUE;
      goto push;
    }
    if (pdata->units+pdata->nonunits==0)
      goto push;
    switch (pdata->op) {
    case 1:
      spheapify(pdata,pdata->vheap[j]);
      elt = pdata->heap[0];
      if (cmp_dltd(pdata->bigf,INTERVAL(elt)))
	goto push;
      return;
    case 2:
      spheapify(pdata,pdata->vheap[j]);
      elt = pdata->heap[0];
      if (cmp_dltd(pdata->bige,INTERVAL(elt)))
	goto push;
      return;
    case 3:
      spheapify(pdata,pdata->vheap[j]);
      elt = pdata->heap[0];
      if (cmp_dltd(pdata->bigf,INTERVAL(elt)) ||
	  cmp_dltd(pdata->bige,INTERVAL(elt)))
	goto push;
      return;
    case 4:
      if (pdata->units+pdata->nonunits <= 1)
	goto push;
      return;
    }
  }
 push:
  fd_enqueue_global(w, global, 0x5/* MINMAX, append*/);   
}

static SP_BOOL
scalar_product_init MAGIC (HIDDEN_PROTO
			   struct linear_data *pdata,
			   SP_BOOL incremental)
{
  int op = pdata->op;
  int i;
  /* Phase 1: compute:
     bigf = rhs - sum{min(a_i x_i)}
     bige = sum{max(a_i x_i)} - rhs

     I_i = max(a_i x_i) - min(a_i x_i)
  */
  if (incremental) {
    for (i=0; i<pdata->ntargets; i++) {
      int j = SV(i);
      Dvar dv = DVAR(j);
      dvar_init(dv, RefAttr(j), RefVar(j));
      if (dvar_is_integer(dv))
	sub_ddd(pdata->rhs,CMIN(j),pdata->rhs);
#if DBG
      {
	long c = COEFF(j);
	struct dword cminj_mem;
	Dword cminj = &cminj_mem;
	struct dword cmaxj_mem;
	Dword cmaxj = &cmaxj_mem;
	TAGGED mint = dvar_min_t(dv);
	TAGGED maxt = dvar_max_t(dv);
	if (c>0) {
	  mul_ltd(c,mint,cminj);
	  mul_ltd(c,maxt,cmaxj);
	  if ((TagIsSmall(mint) && !cmp_deqd(cminj,CMIN(j))) ||
	      (TagIsSmall(maxt) && !cmp_deqd(cmaxj,CMAX(j)))) {
	    SP_ASSERT(FALSE);
	  }
	} else {
	  mul_ltd(c,maxt,cminj);
	  mul_ltd(c,mint,cmaxj);
	  if ((TagIsSmall(maxt) && !cmp_deqd(cminj,CMIN(j))) ||
	      (TagIsSmall(mint) && !cmp_deqd(cmaxj,CMAX(j)))) {
	    SP_ASSERT(FALSE);
	  }
	}
      }
#endif
    }
  } else {
    *pdata->bigf = *pdata->rhs;
    sub_zdd(pdata->bigf,pdata->bige);
    pdata->units = 0;
    pdata->nonunits = 0;
    for (i=0; i<pdata->ntargets; i++) {
      int j = SV(i);
      long c = COEFF(j);
      long ac = ABSCOEFF(j);
      Dvar dv = DVAR(j);
      Dword cminj = CMIN(j);
      Dword cmaxj = CMAX(j);
      Dword intervalj = INTERVAL(j);

      dvar_init(dv, RefAttr(j), RefVar(j));
      if (c>0) {
	mul_ltd(c,dvar_min_t(dv),cminj);
	mul_ltd(c,dvar_max_t(dv),cmaxj);
      } else {
	mul_ltd(c,dvar_max_t(dv),cminj);
	mul_ltd(c,dvar_min_t(dv),cmaxj);
      }
      sub_ddd(cmaxj,cminj,intervalj);
      if (op<4) {
	pdata->heap[i] = j;
	pdata->vheap[j] = i;
      }
      sub_ddd(pdata->bigf,cminj,pdata->bigf);
      add_ddd(pdata->bige,cmaxj,pdata->bige);
      if (dvar_is_integer(dv)) {
	sub_ddd(pdata->rhs,cminj,pdata->rhs);
      } else if (ac==1)
	pdata->units++;
      else
	pdata->nonunits++;
    }
    pdata->heapsize = i;
    if (op<4)
      for (i=(pdata->heapsize-2)>>1; i>=0; i--)
	spheapify(pdata,i);
  }
  return TRUE;
}

static SP_BOOL
scalar_product_le MAGIC (HIDDEN_PROTO
			 struct linear_data *pdata,
			 TERM elt)
{
  long c = COEFF(elt);
  Dvar dv = DVAR(elt);
  struct dword decr_mem;
  Dword decr = &decr_mem;
  struct dword cmax0_mem;
  Dword cmax0 = &cmax0_mem;

  *cmax0 = *CMAX(elt);

  /*
    Phase 2:

    For <=
    ******
    bigf>=0 is a necessary condition.
    bige<=0 is a sufficient condition.

    enforce:
    bigf >= I_i for all i

    rules:
    x_i <=  floor(F / a_i) + min(x_i)  if a_i>0
    x_i >= -floor(F /-a_i) + max(x_i)  if a_i<0

  */

  if (cmp_dltz(pdata->bigf))
    return FALSE;
  if (c>0) {
    struct dword tmp_mem;
    Dword tmp = &tmp_mem;
    struct dword ub_mem;
    Dword ub = &ub_mem;

    add_ddd(pdata->bigf,CMIN(elt),tmp); /* [MC] 3.11.3 */
    div_dldf(tmp,c,ub);
    if (issmall_d(ub)) {/* representable ub */
      if (dvar_fix_max_l(dv,ub->lsw.asword)<0)
	return FALSE;
      mul_ltd(c,dvar_max_t(dv),CMAX(elt));
    } else if (cmp_dgtz(ub)) { /* ub overflow - can't prune the variable */
      long modc = mod_duu(pdata->bigf,c);
      
      add_ddd(pdata->bigf,CMIN(elt),CMAX(elt)); /* ensure bigf >= I_i */
      if (modc>0)
	sub_dud(CMAX(elt),modc,CMAX(elt));
    } else {			/* ub underflow - signal int overflow */
      fd.fd_overflow = TRUE;
      return FALSE;
    }
  } else {
    struct dword tmp_mem;
    Dword tmp = &tmp_mem;
    struct dword lb_mem;
    Dword lb = &lb_mem;

    add_ddd(pdata->bigf,CMIN(elt),tmp); /* [MC] 3.11.3 */
    div_dldc(tmp,c,lb);
    if (issmall_d(lb)) {/* representable lb */
      if (dvar_fix_min_l(dv,lb->lsw.asword)<0)
	return FALSE;
      mul_ltd(c,dvar_min_t(dv),CMAX(elt));
    } else if (cmp_dltz(lb)) {	/* lb underflow - can't prune the variable */
      long modc = mod_duu(pdata->bigf,-c);
      
      add_ddd(pdata->bigf,CMIN(elt),CMAX(elt)); /* ensure bigf >= I_i */
      if (modc>0)
	sub_dud(CMAX(elt),modc,CMAX(elt));
    } else {			/* lb overflow - signal int overflow */
      fd.fd_overflow = TRUE;
      return FALSE;
    }
  }
  sub_ddd(cmax0,CMAX(elt),decr);
  sub_ddd(INTERVAL(elt),decr,INTERVAL(elt));
  sub_ddd(pdata->bige,decr,pdata->bige);
  spheapify(pdata,0);
  if (dvar_is_integer(dv)) {
    sub_ddd(pdata->rhs,CMAX(elt),pdata->rhs);
    if (ABSCOEFF(elt)==1)
      pdata->units--;
    else
      pdata->nonunits--;
  }
  return TRUE;
}

static SP_BOOL
scalar_product_ge MAGIC(HIDDEN_PROTO
			struct linear_data *pdata,
			TERM elt)
{
  long c = COEFF(elt);
  Dvar dv = DVAR(elt);
  struct dword decr_mem;
  Dword decr = &decr_mem;
  struct dword cmin0_mem;
  Dword cmin0 = &cmin0_mem;

  *cmin0 = *CMIN(elt);
  
  /*
    Phase 2:

    For >=
    ******
    bige>=0 is a necessary condition.
    bigf<=0 is a sufficient condition.

    enforce:
    bige >= I_i for all i

    rules:
    x_i >= -floor(E / a_i) + max(x_i)  if a_i>0
    x_i <=  floor(E /-a_i) + min(x_i)  if a_i<0
  */

  if (cmp_dltz(pdata->bige))
    return FALSE;
  if (c>0) {
    struct dword tmp_mem;
    Dword tmp = &tmp_mem;
    struct dword lb_mem;
    Dword lb = &lb_mem;

    sub_ddd(pdata->bige,CMAX(elt),tmp); /* [MC] 3.11.3 */
    div_dldc(tmp,-c,lb);
    if (issmall_d(lb)) {/* representable lb */
      if (dvar_fix_min_l(dv,lb->lsw.asword)<0)
	return FALSE;
      mul_ltd(c,dvar_min_t(dv),CMIN(elt));
    } else if (cmp_dltz(lb)) {	/* lb underflow - can't prune the variable */
      long modc = mod_duu(pdata->bige,c);
      
      sub_ddd(CMAX(elt),pdata->bige,CMIN(elt)); /* ensure bige >= I_i */
      if (modc>0)
	add_dud(CMIN(elt),modc,CMIN(elt));
    } else {			/* lb overflow - signal int overflow */
      fd.fd_overflow = TRUE;
      return FALSE;
    }
  } else {
    struct dword tmp_mem;
    Dword tmp = &tmp_mem;
    struct dword ub_mem;
    Dword ub = &ub_mem;

    sub_ddd(pdata->bige,CMAX(elt),tmp); /* [MC] 3.11.3 */
    div_dldf(tmp,-c,ub);
    if (issmall_d(ub)) {/* representable ub */
      if (dvar_fix_max_l(dv,ub->lsw.asword)<0)
	return FALSE;
      mul_ltd(c,dvar_max_t(dv),CMIN(elt));
    } else if (cmp_dgtz(ub)) {	/* ub overflow - can't prune the variable */
      long modc = mod_duu(pdata->bige,-c);
      
      sub_ddd(CMAX(elt),pdata->bige,CMIN(elt)); /* ensure bige >= I_i */
      if (modc>0)
	add_dud(CMIN(elt),modc,CMIN(elt));
    } else {			/* ub underflow - signal int overflow */
      fd.fd_overflow = TRUE;
      return FALSE;
    }
  }
  sub_ddd(CMIN(elt),cmin0,decr);
  sub_ddd(INTERVAL(elt),decr,INTERVAL(elt));
  sub_ddd(pdata->bigf,decr,pdata->bigf);
  spheapify(pdata,0);
  if (dvar_is_integer(dv)) {
    sub_ddd(pdata->rhs,CMIN(elt),pdata->rhs);
    if (ABSCOEFF(elt)==1)
      pdata->units--;
    else
      pdata->nonunits--;
  }
  return TRUE;
}

/* Preconditions:
   0<coeff<mod, 0=<rhs<mod, 0<mod

   Solve min X such that coeff*X = rhs (modulo mod)
*/
static long solve_gcd(long coeff, long rhs, long mod)
{
  if (rhs==0)
    return 0;
  else {
    long rhs1 = rhs%coeff;
    if (rhs1==0)
      return rhs/coeff;
    else
      return (rhs + mod*solve_gcd(mod%coeff, coeff-rhs1, coeff)) / coeff;
  }
}

/* Preconditions: 
   0<mod, gcdall = gcd(coeff,mod)

   Adjusts minx up and maxx down s.t.
   Returns smallest s.t. coeff*minx = coeff*maxx = rhs (modulo mod)
*/
static void 
adjust_bounds_gcd(long coeff, Dword rhs,
		  long mod, long gcdall,
		  long *minx, long *maxx)
{
  long minx0 = *minx, maxx0 = *maxx;
  long q = mod/gcdall;
  long r, x, s;
  struct dword rhslocal_mem;
  Dword rhslocal = &rhslocal_mem;

  if (coeff>0) {
    *rhslocal = *rhs;
  } else {
    sub_zdd(rhs,rhslocal);
    coeff = -coeff;
  }
  coeff = coeff % mod;
  s = solve_gcd(coeff, (long)mod_duu(rhslocal,mod), mod);
  r = minx0 % q;
  if (r<0)
    r += q;
  x = minx0 - r + s;
  if (x<minx0)
    x += q;
  *minx = x;
  r = maxx0 % q;
  if (r<0)
    r += q;
  x = maxx0 - r + s;
  if (x>maxx0)
    x -= q;
  *maxx = x;
}

SP_BOOL
fd_linear_filter MAGIC(HIDDEN_PROTO Argdecl,SP_term_ref State)
{
  TAGGED tvec, telt, handle, t1;
  SP_BOOL committed; /* TRUE if state can't be backtracked over */
  int nvars, i;
  long state_stamp;
  int ent = -1, total_size;
  int nonground = 0;
  char *ptr;
  struct linear_data *pdata;

  w->numstack_end = NULL;
/*    X(0) = RefTerm(State0); */
  dvar_export_start(w);
  RefTerm(State) = fd_unify_output_state(w,&handle,&state_stamp,&committed);

  if (!IsVar(handle)) {		/* got [Flag | '$free'(Ptr)] */
    pdata = Pdata(struct linear_data,handle);
    nvars = pdata->nvars;
  } else {			/* build persistent state */
    DerefArg(tvec,X(0),1);	/* get CX0 */
    nvars = fd_list_length(tvec);	/* count terms */
    total_size = 
      nvars*sizeof(struct dvar) +
      3*nvars*sizeof(struct dword) +
      3*nvars*sizeof(long) +
      2*nvars*sizeof(TERM) +
      nvars*sizeof(int);
  
    pdata = Palloc(struct linear_data, total_size, handle);
    pdata->rhs = &pdata->rhs_mem;
    pdata->bige = &pdata->bige_mem;
    pdata->bigf = &pdata->bigf_mem;
    ptr = (char *)(pdata+1);
    pdata->term.cmin = (Dword )ptr;
    ptr += nvars*sizeof(struct dword);
    pdata->term.cmax = (Dword )ptr;
    ptr += nvars*sizeof(struct dword);
    pdata->term.interval = (Dword )ptr;
    ptr += nvars*sizeof(struct dword);
    pdata->dvar = (Dvar)ptr;
    ptr += nvars*sizeof(struct dvar);
    pdata->term.coeff = (long *)ptr;
    ptr += nvars*sizeof(long);
    pdata->term.abscoeff = (long *)ptr;
    ptr += nvars*sizeof(long);
    pdata->term.gcd = (long *)ptr;
    ptr += nvars*sizeof(long);
    pdata->target = (TERM *)ptr;
    ptr += nvars*sizeof(TERM);
    pdata->heap = (TERM *)ptr;
    ptr += nvars*sizeof(TERM);
    pdata->vheap = (int *)ptr;
    ptr += nvars*sizeof(int);
    SP_ASSERT(ptr == (char *)(pdata+1)+total_size);
    pdata->overflow = FALSE;
    FD_STORE_SPENV(pdata->spenv);
    pdata->refbase = SP_alloc_globrefs(nvars<<1);
    pdata->stamp = state_stamp-1; /* don't trust initially */
    DerefArg(telt,X(0),2);
    pdata->op = GetSmall(telt) & 0x7;
    pdata->destructor = linear_destructor;
    pdata->daemon = linear_daemon;
    pdata->nvars = nvars;
    pdata->ntargets = nvars;
    DerefArg(tvec,X(0),1);	/* get CX0 */
    DerefArg(telt,X(0),3);	/* get RHS */
    add_ztd(telt,pdata->rhs);
				/* xfer all info to the struct linear_terms */
    for (i=0; i<nvars; i++) {
      TERM elt = i;
      long c;
      
      DerefCar(telt,tvec);
      DerefCdr(tvec,tvec);
      DerefArg(t1,telt,1);
      COEFF(elt) = c = GetSmall(t1);
      ABSCOEFF(elt) = (c>=0 ? c : -c);
      fd_get_var_and_attr(telt+WD(1),RefAttr(elt));
      SV(i) = elt;
    }
    CTagToArg(X(0),1) = atom_nil; /* [MC] 3.12: free for GC */
    for (i=0; i<pdata->ntargets; i++) {
      int j = SV(i);
      Dvar dv = DVAR(j);
      TAGGED functor;
      dvar_init(dv, RefAttr(j), RefVar(j));
      switch (pdata->op) {
      case 4:
	functor = fd.functor_val;
	break;
      default:
	functor = fd.functor_minmax;
	break;
      }
      dvar_attach_daemon(w, dv, pdata, X(1), functor);
    }
  }
  
				/* RESUME HERE */

  fd.fd_overflow = pdata->overflow;
  if (pdata->op==4 || /* #\= case: CMIN/CMAX are not known to be up to date */
      state_stamp != pdata->stamp || /* after backtracking: CMIN/CMAX are not known to be up to date */
      !scalar_product_init(pdata,TRUE)) {
    DerefArg(telt,X(0),3);
    add_ztd(telt,pdata->rhs);
    DerefArg(telt,X(0),4);
    pdata->ntargets = pdata->nvars-GetSmall(telt);
    scalar_product_init(pdata,FALSE);
  }
  pdata->stamp = state_stamp+1;
    
  /* (dis)entailment tests based on GCD */
  switch (pdata->op) {
  case 3:
    refresh_gcd(pdata);	/* compute gcdall, GCD() */
    if (pdata->gcdall>1 && mod_duu(pdata->rhs,pdata->gcdall)!=0)
      return -1;
    break;
  case 4:
    refresh_gcd(pdata);	/* compute gcdall, GCD() */
    if (pdata->gcdall>1 && mod_duu(pdata->rhs,pdata->gcdall)!=0) {
      return 1;
    }
    break;
  }
 loop:
  switch (pdata->op) {
  case 1:			/* #=< */
    if (pdata->heapsize>0) {
      TERM elt = pdata->heap[0];
      if (cmp_dltd(pdata->bigf,INTERVAL(elt))) {
	if (!scalar_product_le(pdata,elt)) {
	  return -1;
	} else {
	  goto loop;
	}
      }
    }
    break;
  case 2:			/* #>= */
    if (pdata->heapsize>0) {
      TERM elt = pdata->heap[0];
      if (cmp_dltd(pdata->bige,INTERVAL(elt))) {
	if (!scalar_product_ge(pdata,elt)) {
	  return -1;
	} else {
	  goto loop;
	}
      }
    }
    break;
  case 3:			/* #= */
    if (pdata->heapsize>0) {
      TERM elt = pdata->heap[0]; 
      if (cmp_dltd(pdata->bigf,INTERVAL(elt))) {
	if (!scalar_product_le(pdata,elt))
	  return -1;
	goto loop;
      }
      if (cmp_dltd(pdata->bige,INTERVAL(elt))) {
	if (!scalar_product_ge(pdata,elt))
	  return -1;
	goto loop;
      }
      refresh_gcd(pdata);	/* compute gcdall, GCD() */
      if (mod_duu(pdata->rhs,pdata->gcdall) != 0)
	return -1;
      if (pdata->units<=1)
	for (i=0; i<pdata->ntargets; i++) {
	  TERM elt = SV(i);
	  Dvar dv = DVAR(elt);
	  TAGGED tmin = dvar_min_t(dv);
	  TAGGED tmax = dvar_max_t(dv);

	  if (AreSmall(tmin,tmax) && tmin!=tmax && GCD(elt)>pdata->gcdall) {
	    /* Ensure that:

	       Ai * min(Xi) - RHS = Ai * max(Xi) - RHS = 0 (modulo G)

	       where G is the gcd of all coefficients excepts Ai.
	       Fails if an empty domain is produced.
	    */
	    long imin = GetSmall(tmin);
	    long imax = GetSmall(tmax);
	    int rc;

	    adjust_bounds_gcd(COEFF(elt),
			      pdata->rhs,
			      GCD(elt),
			      pdata->gcdall,
			      &imin, &imax);
	    rc = dvar_fix_interval_l(dv, imin, imax);
	    if (rc<0)
	      return -1;
	    if (rc>0) {
	      struct dword cmin0_mem;
	      Dword cmin0 = &cmin0_mem;
	      struct dword cmax0_mem;
	      Dword cmax0 = &cmax0_mem;
	      struct dword decr_mem;
	      Dword decr = &decr_mem;
	      long c = COEFF(elt);

	      *cmin0 = *CMIN(elt);
	      *cmax0 = *CMAX(elt);

	      if (c>0) {
		mul_ltd(c,dvar_min_t(dv),CMIN(elt));
		mul_ltd(c,dvar_max_t(dv),CMAX(elt));
	      } else {
		mul_ltd(c,dvar_max_t(dv),CMIN(elt));
		mul_ltd(c,dvar_min_t(dv),CMAX(elt));
	      }
	      sub_ddd(CMAX(elt),CMIN(elt),INTERVAL(elt));
	      sub_ddd(cmax0,CMAX(elt),decr);
	      sub_ddd(pdata->bige,decr,pdata->bige);
	      sub_ddd(CMIN(elt),cmin0,decr);
	      sub_ddd(pdata->bigf,decr,pdata->bigf);
	      spheapify(pdata,pdata->vheap[elt]);
	      if (dvar_is_integer(dv)) {
		sub_ddd(pdata->rhs,CMIN(elt),pdata->rhs);
		if (ABSCOEFF(elt)==1)
		  pdata->units--;
		else
		  pdata->nonunits--;
	      }
	      goto loop;
	    }
	  }
	}
    }
    break;
  case 4:			/* #\= */
    nonground = pdata->units+pdata->nonunits;
    if (nonground==0 && cmp_deqz(pdata->rhs))
      return -1;
    else if (nonground==1) {
      for (i=0; cmp_deqz(INTERVAL(SV(i))); i++)
	;
      {
	TERM elt = SV(i);
	Dvar dv = DVAR(elt);

	if (mod_duu(pdata->rhs,ABSCOEFF(elt))==0) { /* RHS a multiple of coefficient */
	  dvar_prune_value_l(dv,div_dll(pdata->rhs,COEFF(elt)),w);
	}
      }
    }
    break;
  }
  switch (pdata->op) {
  case 1:
    ent = (cmp_dltz(pdata->bigf) ? -1 : cmp_dlez(pdata->bige));
    break;
  case 2:
    ent = (cmp_dltz(pdata->bige) ? -1 : cmp_dlez(pdata->bigf));
    break;
  case 3:
    ent = (cmp_dltz(pdata->bigf) ? -1 :
	   cmp_dltz(pdata->bige) ? -1 : 
	   cmp_dgtz(pdata->bigf) ? 0 :
	   cmp_dgtz(pdata->bige) ? 0 : 1);
    break;
  case 4:
    ent = (nonground<=1);
    break;
  }

  if (ent>=0) {
    int inf = 0;
    int sup = pdata->ntargets-1;
    TERM held = SV(sup); /* sup is the hole */
    TERM current = SV(inf);
    
    while (inf<=sup) {
      Dvar dv = DVAR(current);
      
      dvar_pruning_done(w,dv);
      dvar_export(w,dv);
      if (!dvar_is_integer(dv)) {
	SV(inf) = current;
	inf++;
	current = (inf>=sup ? held : SV(inf));
      } else {
	SV(sup) = current;
	sup--;
	current = (inf>=sup ? held : SV(sup));
      }
    }
    pdata->ntargets = inf;
  }
  if (ent==0 && pdata->op==3 && pdata->ntargets==2 && cmp_deqz(pdata->rhs)) {
    TERM elt0 = SV(0);
    TERM elt1 = SV(1);

    if (COEFF(elt0) == -COEFF(elt1)) {
      dvar_export_equal(w,DVAR(elt0), DVAR(elt1));
      ent = 1;
    }
  }
  if (islong_d(pdata->rhs)) {
    CTagToArg(X(0),3) = MakeInteger(pdata->rhs->lsw.asword);
  } else if (isfinite_d(pdata->rhs)) {			/* make a two-word bignum */
    TAGGED *h;
    RequireHeap(4,EVAL_ARITY);
    h = w->global_top;
    h[0] = BIGNUM_HEADER-LStep(1);
    h[1] = pdata->rhs->lsw.asuword;
    h[2] = pdata->rhs->msw.asuword;
    h[3] = BIGNUM_HEADER-LStep(1);
    w->global_top = h+4;
    CTagToArg(X(0),3) = MakeStructure(h);
  } else {			/* RHS cannot be represented in two-word bignum */
    fd.fd_overflow = TRUE;
    ent = -1;
  }
  CTagToArg(X(0),4) = MakeSmall(pdata->nvars-pdata->ntargets);
  return ent;
}
